function ROI1Profile = ShearView(SV_parameters)

% FIXME: IMPORTANT: start program in patient root directory, SHEARRATE WILL WRITE
% FILES AND FOLDERS IN THIS DIRECTORY 
% written by Alex Barker (with code modified from Craig Lanning, TCH)

% pull out file info from SVparameters structure
IntDirName  = SV_parameters.Intdir;
IntFileName = SV_parameters.Intfile;
PCDirName   = SV_parameters.Veldir;
PCFileName  = SV_parameters.Velfile;
ARGUSdir    = SV_parameters.ROIdir;
ARGUSfile   = SV_parameters.ROIfile;
WkgDir      = SV_parameters.wkgdir;
WkspaceName = SV_parameters.wkspaceName;
PatientID   = SV_parameters.PatientID;
SystoleMeth = SV_parameters.Systole;
Platform    = SV_parameters.Platform;
if isfield(SV_parameters,'DBname')
    DBtype  = SV_parameters.DBname;
else
    DBtype  = 'test';
end
intrpmeth   = SV_parameters.interp(2:end);

% What should be plotted and saved?      
VelFlowFLAG     = SV_parameters.saveVelFlow;    % Save Velocity Profile to tiff file
movieFLAG       = SV_parameters.movVelFlow;     % write movie file
PlotAllROIs     = 0;                            % need this for Get ROI's (may change this...)

% set working directory as the active path
cd(WkgDir);

% Determine if we are dealing with a Segment or Argus file:
[pathstr name ext] = fileparts([ARGUSdir '\' ARGUSfile]); clear pathstr name versn
if strcmp(ext,'.txt')
    ARGUS_ROI = 1;
elseif strcmp(ext,'.con') % this is a medis ROI
    ARGUS_ROI = 3;
elseif strcmp(ext,'.seg') % this is a segment ROI
    ARGUS_ROI = 4;
else
    display('error in ROI file type extension: .txt, .seg, or .con file not detected')
    display('    *.txt - Segment contour')
    display('    *.seg - Segment contour')
    display('    *.con - Medis contour')
    return
end

ZoomOut             = 15;           % Zoom factor - out from ROI
Elevation           = 40.0;         % Elevation for 3D plots

scrsz = get(0,'ScreenSize');        % Determine the Overall Screen Size and 
scrsz(4) = scrsz(4)-50;             % adjust for windows taskbar
% *************************************************************************

% Check File Agreement 
[countInt,SeriesNumInt] = countDICOMFiles(IntDirName);
[countPC,SeriesNumPC]   = countDICOMFiles(PCDirName);

if countInt ~= countPC
    
    disp('Error: Incorrect number of DICOM images.');
%    return;
    disp('...but continuing anyway, in case of 3dir sequence');
    
end

% First Intensity Info
filepath                = [IntDirName '\' IntFileName];  % This is assuming the path is already set
IntensityInfo           = dicominfo(filepath);
structInt               = getDICOMinfo(IntensityInfo);
TempIntensityImageData  = dicomread(IntensityInfo);


% First Velocity component info.
filepath                = [PCDirName '\' PCFileName];   
VelDirOneInfo           = dicominfo(filepath);              
structPC                = getDICOMinfo(VelDirOneInfo);       
TempVelDirOneImageData  = dicomread(VelDirOneInfo);

MagnetTT(1)            = VelDirOneInfo.TriggerTime; % note the TT is for PC image!!

% Velocity Encoding - get VENC value

if Platform == 1 % Siemens
    % Find the '_v' tag for the beginning (for Siemens Dicoms)
    venc_startindex = strfind(VelDirOneInfo.SequenceName,'_v');
    venc_endindex   = strfind(VelDirOneInfo.SequenceName,'in');

    % Check for inplane directions
    if isempty(venc_endindex)
        disp('WARNING! Venc may not be correct.')
    end
    VENC = str2num(VelDirOneInfo.SequenceName(venc_startindex+2:venc_endindex-1));
    clear('venc_startindex','venc_endindex');
elseif Platform == 2 % Philips
    VENC = 100; %TODO: need to fix this!!
elseif Platform == 3 % GE Healthcare
    VENC = 200; %TODO: need to fix!!
end

if Platform ==3
    TempVelDirOneVelData = (double(TempVelDirOneImageData))*VENC/2048;   
else
    TempVelDirOneVelData = (double(TempVelDirOneImageData)-2048)*VENC/2048;  
end

%--------ROI Calculation (performed on first intensity image)-----

% Determine whether to use ARGUS ROI, Medis or Segment to complete
% flow analysis

if ARGUS_ROI == 1
   
    % ARGUS generated ROI coodinates
    [ROI t numReg ROInumber ROIfigure] = GetARGUSrois(TempIntensityImageData(:,:,1), [ARGUSdir '/' ARGUSfile], PlotAllROIs);
    
    for i = 1:t
        
        % pull ROI throughout cine loop
        ROI_X = ROI(:,1,i,ROInumber);  
        ROI_Y = ROI(:,2,i,ROInumber);  
        b(:,:,i) = [ROI_X ROI_Y];
        BW(:,:,i) = poly2mask(ROI_X, ROI_Y, ...  %% Create to calc xmin, ymin...
            size(TempIntensityImageData,1),...
            size(TempIntensityImageData,2));
    
    end
    
elseif ARGUS_ROI == 3
    medisContour = MedisContourLoad([ARGUSdir '\' ARGUSfile]);
    if medisContour.Zoom ~= 1
        disp('Warning: Contour scaling may not be correct!');
    end
    ROInumber = 1;
    t = medisContour.Curve(end).Phase+1;
    for i = 1:t
        temp_contour(:,:) = [medisContour.Curve(i).Contour_U medisContour.Curve(i).Contour_V]; % problem is that b changes in size from time pt to time pt... need to interpret to ~200
        temp_length = length(temp_contour(:,1));
        ROI_X     = interp1(1:temp_length,temp_contour(:,1),linspace(1,temp_length,200),'pchip')'; % interpret to 200 pts
        ROI_Y     = interp1(1:temp_length,temp_contour(:,2),linspace(1,temp_length,200),'pchip')'; % interpret to 200 pts
        b(:,:,i) = [ROI_X ROI_Y];
        BW(:,:,i) = poly2mask(ROI_X, ROI_Y, ...  %% Create to calc xmin, ymin...
            size(TempIntensityImageData,1),...
            size(TempIntensityImageData,2));
        clear temp_contour temp_length 
    end
    
    
elseif ARGUS_ROI == 4
   
    % Segment generated ROI coodinates
    [ROI_Y ROI_X] = getSegmentROI([ARGUSdir '\' ARGUSfile]);
    
    ROInumber = 1;
    t = size(ROI_X,2);
    for i = 1:t
        
        % pull ROI throughout cine loop
        b(:,:,i) = [ROI_X(:,i) ROI_Y(:,i)];
        BW(:,:,i) = poly2mask(ROI_X(:,i), ROI_Y(:,i), ...  %% Create to calc xmin, ymin...
            size(TempIntensityImageData,1),...
            size(TempIntensityImageData,2));
    
    end
  
end

% The following if statement looks for temporal ROI by look at ARGUS_ROI, if it is a
% temporally variable ROI, then Xmin, max.... are calculated differently

[Xmin Xmax Ymin Ymax IntensityImageData VelDirOneVelData BW L stats b CropRegion XTextSpacer YTextSpacer]...
        = CropAndProp(TempIntensityImageData, TempVelDirOneVelData, BW, countInt, ZoomOut, t, b);
    
% LOAD DATA --------------------------------------------------------------%

IntImStart = str2double(IntFileName(9:12));
PCImStart  = str2double(PCFileName(9:12));
[IntensityImageData VelDirOneVelData MagnetTT] = LoadSequence(countInt, ...
    IntDirName, SeriesNumInt, IntImStart, PCDirName, SeriesNumPC, PCImStart,...
    IntensityImageData, VelDirOneVelData, CropRegion, VENC, MagnetTT, ROInumber,Platform);

% FIND STATS -------------------------------------------------------------%
% calculate uncertainty at each pixel ------------------------------------%

[VelDironeVelDataMask VelDironeVelDataAverage VelDironeVelDataUncertainty...
    flowRateVsTime averageUncertainty averageVelocity avgVelocityVsTime ...
    peakVelocityVsTime CrossSectionArea FlowRate]...
    = FlowStats(VelDirOneVelData,ARGUS_ROI, BW, structInt.PixelSpacingRows,...
    structInt.PixelSpacingColumns, t, ROInumber);
        
% correct for negative velocity values
tempfigure = figure;
plot(flowRateVsTime)
DirChangeReply = questdlg('Do you want to change the direction of the Flow?','Option to Change Flow','No');
if strcmp(DirChangeReply,'Yes')
    VelDironeVelDataMask = -VelDironeVelDataMask;
    flowRateVsTime = -flowRateVsTime;
    averageVelocity = -averageVelocity;
    avgVelocityVsTime = -avgVelocityVsTime;
    FlowRate = -FlowRate;
end
close(tempfigure)
    
    
% find max flow rate (flow systole) & max/min area
MaxFlow = max(flowRateVsTime);

tempcell = struct2cell(stats);
tempmat  = cell2mat(tempcell(1,:,:));
Areas   = tempcell(1,:);
Areas   = cell2mat(Areas);
Areas   = Areas*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*.01;
MaxArea = max(tempmat);
MaxArea = MaxArea*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*.01;
MinArea = min(tempmat);
MinArea = MinArea*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*.01;

[positiveFlow negativeFlow] = flowRatio(VelDironeVelDataMask);
positiveFlow = positiveFlow*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*0.0006;
negativeFlow = negativeFlow*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*0.0006;
if SystoleMeth
    [MaxFlowMagnetud SystoleFrame] = max(abs(flowRateVsTime));
else
    [MaxPositiveFlow SystoleFrame] = max(positiveFlow);
end

% Find Maximum/minimum velocity and velocity at systole
VelMax           = max(VelDironeVelDataMask(:,:,:));
VelMax           = max(max(VelMax));
[VelMaxSys Cmax] = max(VelDironeVelDataMask(:,:,SystoleFrame));
[VelMaxSys Rmax] = max(VelMaxSys);

VelMin           = min(VelDironeVelDataMask(:,:,:));
VelMin           = min(min(VelMin));
[VelMinSys Cmin] = min(VelDironeVelDataMask(:,:,SystoleFrame));
[VelMinSys Rmin] = min(VelMinSys);




%% Flow/Function/Profile Comparison through Centroid

% Generate the initial figure
fig_width  = 800;
fig_height = 400;
count = countInt;

if VelFlowFLAG || movieFLAG
    ROI1Profile = figure('Position', [scrsz(3)-835 scrsz(4)-fig_height 675 200*2],...
        'Color',[1 1 1], 'Name',['ROI ' num2str(ROInumber) ' Profile Comparison through Centroid' ],...
        'toolbar','none','Menubar','none');

    aleft = 0;
    abottom = 0;
    awidth = .34;
    aheight = .5;
    aindex = 0.33;

    ax1  = axes;
    set(ax1,'OuterPosition', [aleft abottom+aheight awidth aheight])

    %PLOT initial image of intensity map-------------
    % find indexes for nonnegative values of the ROI boundary coords
    y = find(b(:,1,1) > 0);
    x = find(b(:,2,1) > 0);

    mapdepth = 256;
    vmap = cat(1, jet(mapdepth/2), gray(mapdepth/2));
    colormap(vmap);

    hInt          = image(IntensityImageData(:,:,1),'CDataMapping', 'scaled');
    axis image; hold on;
    hROI          = plot(b(y,1,1),b(x,2,1),'g','LineWidth',2);
    hCentroid     = plot(stats(1).Centroid(1), stats(1).Centroid(2), '+'); hold off
    htxtXCentroid = text((Xmin+XTextSpacer),Ymax-YTextSpacer,['Centroid x   = '...
        num2str(stats(1).Centroid(1),3) ' pix'], 'Color', [1 1 1], 'FontSize',8);
    htxtYCentroid = text((Xmin+XTextSpacer),Ymax-2*YTextSpacer,['Centroid y   = '...
        num2str(stats(1).Centroid(2),3) ' pix'], 'Color', [1 1 1], 'FontSize',8);
    htxtArea      = text((Xmin+XTextSpacer),Ymin+2*YTextSpacer,['Area = '...
        num2str(stats(1).Area*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*.01,3) ' cm^2'], ...
        'Color', [1 1 1], 'FontSize',8);
    title('Intensity Image','FontSize',10);
    ylabel('Pixels','FontSize',10);
    xlabel('Pixels','FontSize',10);
    text((Xmin+XTextSpacer),Ymin+1*YTextSpacer,['Region Chosen = '...
        num2str(ROInumber)], 'Color', [1 1 1], 'FontSize',8);
    hold off
    %---------------------------------------------------

    ax2  = axes;
    set(ax2,'OuterPosition', [aleft abottom awidth aheight])

    %PLOT initial image of Velocity map-------------
    hVel = image(VelDironeVelDataMask(:,:,1),'CDataMapping', 'scaled');
    axis image; hold on

    htxtVavg  = text(Xmin+XTextSpacer,Ymax-1*YTextSpacer,['Avg. Vel. = ' num2str(avgVelocityVsTime(1),3) ' cm/s'], 'Color', [1 1 1], 'FontSize',8);
    htxtVpeak = text(Xmin+XTextSpacer,Ymax-2*YTextSpacer,['Peak Vel. = ' num2str(peakVelocityVsTime(1),3) ' cm/s'], 'Color', [1 1 1], 'FontSize',8);
    htxtQ     = text(Xmin+XTextSpacer,Ymin+YTextSpacer,['Q = ' num2str(flowRateVsTime(1),3) ' L/min'], 'Color', [1 1 1], 'FontSize',8);
    text(Xmin+XTextSpacer,Ymin+2*YTextSpacer,['Qavg = ' num2str(FlowRate,3) ' L/min'], 'Color', [1 1 1], 'FontSize',8);
    title('Velocity Map','FontSize',10);
    ylabel('Pixels'     ,'FontSize',10);
    xlabel('Pixels'     ,'FontSize',10);

    hlineHor  = line([Xmin Xmax], [stats(1).Centroid(2) stats(1).Centroid(2)],'linewidth',2,'Linestyle','--','Color','b');
    hlineVert = line([stats(1).Centroid(1) stats(1).Centroid(1)], [Ymin Ymax],'linewidth',2,'Linestyle','--','Color','r');
    hold off;

    CmLength   = size(get(gcf,'Colormap'),1);% Colormap length
    CLim1      = get(ax1,'CLim');            % CLim values for each axis
    set(ax1,'CLim',newclim(mapdepth/2,mapdepth,CLim1(1)-2,CLim1(2),CmLength))
    set(ax2,'CLim',newclim(1,mapdepth,-VENC/1.5,VENC/1.5,CmLength))
    colorbar('OuterPosition',[.279 .048 .071 .423])

    %---------------------------------------------------

    haxes = axes;
    set(haxes,'OuterPosition', [aleft+1.05*aindex abottom awidth aheight])

    %PLOT initial Velocity Profile-----------------------
    [m n] = size(VelDironeVelDataMask(:,:,1));
    UnrefinedX = (1:n);
    UnrefinedY = (1:m);
    UnrefinedVelProfileX(1,:) = VelDironeVelDataMask(round(stats(1).Centroid(2)),:,1);
    UnrefinedVelProfileY(:,1) = VelDironeVelDataMask(:,round(stats(1).Centroid(1)),1);
    hlineXprof = plot(UnrefinedX, UnrefinedVelProfileX(1,:),'LineWidth',1); hold on;
    hlineYprof = plot(UnrefinedY, UnrefinedVelProfileY(:,1),'r','LineWidth',1); hold off;
    if (Xmax-Xmin) >= (Ymax-Ymin)
        axis([Xmin Xmax -VENC VENC]);
        htxtTT   = text(Xmin+XTextSpacer,VENC-VENC*.2,['Trigger time= ' num2str(MagnetTT(1),4) ' ms'], 'Color', [0 0 0], 'FontSize',8);
        htxtVenc = text(Xmin+XTextSpacer,-VENC+VENC*.2,['VENC= ' num2str(VENC,3) ' ms'], 'Color', [0 0 0], 'FontSize',8);
    else
        axis([Ymin Ymax -VENC VENC]);
        htxtTT   = text(Ymin+YTextSpacer,VENC-VENC*.2,['Trigger time= ' num2str(MagnetTT(1),4) ' ms'], 'Color', [0 0 0], 'FontSize',8);
        htxtVenc = text(Ymin+YTextSpacer,-VENC+VENC*.2,['VENC= ' num2str(VENC,3) ' ms'], 'Color', [0 0 0], 'FontSize',8);
    end
    title('Velocity (along dashed lines)','FontSize',10);
    ylabel('Velocity [cm/s]','FontSize',10);
    xlabel('Pixels','FontSize',10);
    %---------------------------------------------------

    flowaxes = axes;
    set(flowaxes, 'OuterPosition', [aleft+.9*aindex abottom+aheight awidth*2.1 aheight])

    %PLOT initial Flow Graph      -----------------------

    plot(MagnetTT/1000,negativeFlow,'Color','red', 'LineWidth',1);
    hold on
    plot(MagnetTT/1000,positiveFlow,'Color','g','LineWidth',1);
    plot(MagnetTT/1000,flowRateVsTime,'Color','b','Marker','o','LineWidth',2);
    title( 'Flow Rate vs Time','FontSize',10);
    xlabel('Time [s]'         ,'FontSize',10);
    ylabel('Flow Rate [L/min]','FontSize',10);
    grid on
    axis([0 max(MagnetTT)/1000 min(negativeFlow)*1.05 max(positiveFlow)*1.05])
    hlinetime = plot([MagnetTT(1)/1000 MagnetTT(1)/1000],[min(negativeFlow)*1.05 max(positiveFlow)*1.05],'Color','r','LineWidth',2);
    hold off
    leg = legend('Negative Flow', 'Positive Flow', 'Net Flow Rate');
    set(leg,'EdgeColor',[1 1 1],'FontSize',10)

    %---------------------------------------------------
    defaxes = axes;
    set(defaxes,'OuterPosition', [aleft+2*aindex abottom awidth aheight])

    plot(0,0,'Color',[.5 .5 .5],'Marker','+')
    axis([Xmin-stats(1).Centroid(1) Xmax-stats(1).Centroid(1) Ymin-stats(i).Centroid(2)+4 Ymax-stats(i).Centroid(2)+4]);
    axis ij; axis equal; axis manual; box on; hold on
    for i = 1:count
        % find indexes for nonnegative values of the ROI boundary coords
        y = find(b(:,1,i) > 0);
        x = find(b(:,2,i) > 0);
        plot(b(y,1,i)-stats(i).Centroid(1) ,b(x,2,i)-stats(i).Centroid(2) ,'Color',[.5 .5 .5],'LineWidth',2);
    end
    hlineDef = plot(b(y,1,1)-stats(1).Centroid(1) ,b(x,2,1)-stats(1).Centroid(2) ,'Color','k','LineWidth',2);
    h = findobj(gca,'Type','line','-and','-not',{'Marker','+'});
    legend(...
        h(1:2),{'Current Position','Deformation'},...
        'FontSize',8,...
        'Position',[0.7472 0.098 0.1985 0.095],...
        'XColor',[1 1 1],...
        'YColor',[1 1 1],...
        'Color','none',...
        'Box', 'Off');
    title( 'Wall Deformation' ,'FontSize',10);
    ylabel('Position [Pixels]','FontSize',10);
    xlabel('Position [Pixels]','FontSize',10);

    PlotCompass;

end

if movieFLAG;
    % save intial frame of movie
    if strcmp(DBtype,'Phantom')
        movfile = '';
    else
        movfile = strrep(WkspaceName, '.mat', '');    
    end
    mkdir('_movies')
    newpath = pwd;
    cd([newpath '\_movies']);
%     mov = avifile([PatientID '_' movfile '_VelFlowProf.avi'],'compression','Cinepak');
    mov = VideoWriter([PatientID '_' movfile '_VelFlowProf.avi'],'Motion JPEG AVI');
    open(mov)

    tempframe = getframe(ROI1Profile);
%     mov = addframe(mov,tempframe);
    writeVideo(mov,tempframe);
    
    for i = 2:1:count

        % find indexes for nonnegative values of the ROI boundary coords
        y = find(b(:,1,i) > 0);
        x = find(b(:,2,i) > 0);

        set(ROI1Profile,'CurrentAxes',ax1);
        set(hInt,'CData',IntensityImageData(:,:,i))
        set(hROI,'XData',b(y,1,i),'YData',b(x,2,i));
        set(hCentroid,'XData',stats(i).Centroid(1),'YData',stats(i).Centroid(2))
        set(htxtXCentroid,'String',['Centroid x   = '...
            num2str(stats(i).Centroid(1),3) ' pix'])
        set(htxtYCentroid,'String',['Centroid y   = '...
            num2str(stats(i).Centroid(2),3) ' pix'])
        set(htxtArea,'String',['Area = '...
            num2str(stats(i).Area*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*.01,3) ' cm^2'])

        set(ROI1Profile,'CurrentAxes',ax2);
        set(hVel,'CData',VelDironeVelDataMask(:,:,i))
        set(htxtVavg ,'String',['Avg. Vel. = ' num2str(avgVelocityVsTime(i),3) ' cm/s'])
        set(htxtVpeak,'String',['Peak Vel. = ' num2str(peakVelocityVsTime(i),3) ' cm/s'])
        set(htxtQ    ,'String',['Q = ' num2str(flowRateVsTime(i),3) ' L/min'])
        set(hlineHor,'XData',[Xmin Xmax],'YData',[stats(i).Centroid(2) stats(i).Centroid(2)])
        set(hlineVert,'XData',[stats(i).Centroid(1) stats(i).Centroid(1)],'YData',[Ymin Ymax])

        % Plot Velocity Profile
        set(ROI1Profile,'CurrentAxes',haxes);
%         [m n] = size(VelDironeVelDataMask(:,:,i));
%         UnrefinedX = [1:n];
%         UnrefinedY = [1:m];
        UnrefinedVelProfileX(i,:) = VelDironeVelDataMask(round(stats(i).Centroid(2)),:,i);
        UnrefinedVelProfileY(:,i) = VelDironeVelDataMask(:,round(stats(i).Centroid(1)),i);
        set(hlineXprof,'XData',UnrefinedX,'YData',UnrefinedVelProfileX(i,:))
        set(hlineYprof,'XData',UnrefinedY,'YData',UnrefinedVelProfileY(:,i))
        set(htxtTT  ,'String',['Trigger time= ' num2str(MagnetTT(i),4) ' ms'])
        set(htxtVenc,'String',['VENC= ' num2str(VENC,3) ' ms'])

        set(ROI1Profile,'CurrentAxes',flowaxes);
        set(hlinetime,'XData',[MagnetTT(i)/1000 MagnetTT(i)/1000],'YData',[min(negativeFlow)*1.05 max(positiveFlow)*1.05])

        set(ROI1Profile,'CurrentAxes',defaxes);
        set(hlineDef,'XData',b(y,1,i)-stats(i).Centroid(1),'YData',b(x,2,i)-stats(i).Centroid(2))

        % prep for movie
        tempframe = getframe(ROI1Profile);
        writeVideo(mov,tempframe);

    end


    close(mov);
    clear newpath mov CI;
    cd('..')

end


%% Plot ROI1 Function at Systole

% Set i = systoleframe
if VelFlowFLAG
    
    i = SystoleFrame;
    y = find(b(:,1,i) > 0);
    x = find(b(:,2,i) > 0);

    set(ROI1Profile,'CurrentAxes',ax1);
    set(hInt,'CData',IntensityImageData(:,:,i))
    set(hROI,'XData',b(y,1,i),'YData',b(x,2,i));
    set(hCentroid,'XData',stats(i).Centroid(1),'YData',stats(i).Centroid(2))
    set(htxtXCentroid,'String',['Centroid x   = '...
        num2str(stats(i).Centroid(1),3) ' pix'])
    set(htxtYCentroid,'String',['Centroid y   = '...
        num2str(stats(i).Centroid(2),3) ' pix'])
    set(htxtArea,'String',['Area = '...
        num2str(stats(i).Area*structInt.PixelSpacingRows*structInt.PixelSpacingColumns*.01,3) ' cm^2'])

    set(ROI1Profile,'CurrentAxes',ax2);
    set(hVel,'CData',VelDironeVelDataMask(:,:,i))
    set(htxtVavg ,'String',['Avg. Vel. = ' num2str(avgVelocityVsTime(i),3) ' cm/s'])
    set(htxtVpeak,'String',['Peak Vel. = ' num2str(peakVelocityVsTime(i),3) ' cm/s'])
    set(htxtQ    ,'String',['Q = ' num2str(flowRateVsTime(i),3) ' L/min'])
    set(hlineHor,'XData',[Xmin Xmax],'YData',[stats(i).Centroid(2) stats(i).Centroid(2)])
    set(hlineVert,'XData',[stats(i).Centroid(1) stats(i).Centroid(1)],'YData',[Ymin Ymax])

    % Plot Velocity Profile
    set(ROI1Profile,'CurrentAxes',haxes);
%     [m n] = size(VelDironeVelDataMask(:,:,i));
%     UnrefinedX = [1:n];
%     UnrefinedY = [1:m];
    UnrefinedVelProfileX(i,:) = VelDironeVelDataMask(round(stats(i).Centroid(2)),:,i);
    UnrefinedVelProfileY(:,i) = VelDironeVelDataMask(:,round(stats(i).Centroid(1)),i);
    set(hlineXprof,'XData',UnrefinedX,'YData',UnrefinedVelProfileX(i,:))
    set(hlineYprof,'XData',UnrefinedY,'YData',UnrefinedVelProfileY(:,i))
    set(htxtTT  ,'String',['Trigger time= ' num2str(MagnetTT(i),4) ' ms'])
    set(htxtVenc,'String',['VENC= ' num2str(VENC,3) ' ms'])

    set(ROI1Profile,'CurrentAxes',flowaxes);
    set(hlinetime,'XData',[MagnetTT(i)/1000 MagnetTT(i)/1000],'YData',[min(negativeFlow)*1.05 max(positiveFlow)*1.05])

    set(ROI1Profile,'CurrentAxes',defaxes);
    set(hlineDef,'XData',b(y,1,i)-stats(i).Centroid(1),'YData',b(x,2,i)-stats(i).Centroid(2))
    

end

%% Save to .mat files ------------------------------------%

% save workspace variables for ROI#1
mkdir('_wkspce')
save(['_wkspce\' intrpmeth '_' WkspaceName]);
display(['Matlab workspace saved as ' pwd '\_wkspace\' WkspaceName]);

% TODO: need to fix figure handles so that they are more robust... ie changing
% the order of the figrues won't affect the print outs

if VelFlowFLAG
    Plotfile = strrep(WkspaceName, '.mat', '');
    mkdir('_flowCalcs')
    set(ROI1Profile,'PaperPositionMode','auto');
    print(ROI1Profile,'-zbuffer', '-dtiff', ['_flowCalcs\' PatientID '_' Plotfile '_VelFlowProf.tif']);
    display(['Velocity & Flow Figure saved as ' pwd '\_flowCalcs\' PatientID '_' Plotfile '_VelFlowProf.tif']);
else
    ROI1Profile = [];
end





